##replication of figure 1##



##replication of panel a##
a=10000;
b=55500;
zstar=10000;
x0=500;
binwidth=x0;

y0=5;

binned_data03070<-binned_data03070[-1,];  ##get rid of the first bin with no loan obs; similar results if the first bin is kept##

bins_excl_r <- 1
         
            firstpass_prep <- bunching::prep_data_for_fit(binned_data03070,zstar =zstar,  binwidth = x0,
                          bins_l =a/x0, bins_r=b/x0,rn=c(5000,10000), bins_excl_r=bins_excl_r,bins_excl_l=0/x0,poly=y0)
            firstpass <- bunching::fit_bunching(firstpass_prep$data_binned, 
                firstpass_prep$model_formula, binwidth=x0, notch=TRUE)
            B_below <- firstpass$B_zl_zstar
            M_above <- -firstpass$B_zstar_zu
            if (M_above > B_below) {
                stop("Missing mass above zstar is larger than bunching mass below. Are you sure this is a notch?")
            }
            available_bins <- (max(firstpass_prep$data_binned$bin) - 
                zstar)/binwidth
            DoF_remaining <- nrow(firstpass_prep$data_binned) - 
                nrow(firstpass$coefficients) - 1
            notch_iterations_bound <- min(available_bins, DoF_remaining)
            zu_bin <- bins_excl_r
            tmp_firstpass_prep <- firstpass_prep
            while ((B_below > M_above) & (zu_bin < notch_iterations_bound)) {
                zu_bin_final <- zu_bin
                zu_bin <- zu_bin + 1
                newvar <- paste0("bin_excl_r_", zu_bin)
                tmp_firstpass_prep$data_binned[[newvar]] <- ifelse(tmp_firstpass_prep$data_binned$z_rel == 
                  zu_bin, 1, 0)
                   tmp_firstpass_prep$model_formula <- stats::as.formula(paste(Reduce(paste, 
                  deparse(tmp_firstpass_prep$model_formula)), 
                  newvar, sep = " + "))
                tmp_firstpass <- bunching::fit_bunching(tmp_firstpass_prep$data_binned, 
                  tmp_firstpass_prep$model_formula, binwidth=x0, notch=TRUE)
                B_below <- tmp_firstpass$B_zl_zstar
                M_above <- -tmp_firstpass$B_zstar_zu
                if (B_below > M_above) {
                  firstpass_prep <- tmp_firstpass_prep
                  firstpass <- tmp_firstpass
                  zu_bin_final <- zu_bin
                }
            }

mb<-zu_bin*x0+zstar;  

plot_binned_data0307<-tmp_firstpass_prep$data_binned;

plot_binned_data0307$cf<-tmp_firstpass$cf_density;

plot_binned_data0307$pfreq<-100*plot_binned_data0307$freq/sum(plot_binned_data0307$freq);
plot_binned_data0307$pcf<-100*plot_binned_data0307$cf/sum(plot_binned_data0307$freq);


pdf(width = 8, # The width of the plot in inches
    height = 4) # The height of the plot in inches


ggplot(plot_binned_data0307)+geom_line(aes(x=bin,y=pfreq),color="black",size=0.4)+geom_line(aes(x=bin,y=pcf),linetype="longdash",color="red",size=0.4)+ 
  labs(x="Loan Amount",y="Percentage of Loans")+geom_vline(xintercept = 9900, linetype="dotted", color = "black", size=0.3)+
scale_x_continuous(breaks=c(5000,15000,25000,35000,45000,55000,65000))+
geom_text(parse=TRUE,size=3,aes(x=9300, label ="italic(K[L])", y=6.5), colour="black", angle=0, text=element_text(size=0.3))+
geom_vline(xintercept = mb, linetype="dotted", 
                color = "black", size=0.3)+
geom_text(parse=TRUE,size=3,aes(x=16700, label ="italic(K[U])" , y=6.5), colour="black", angle=0, text=element_text(size=0.3))+ 
theme_bw() +
  theme(axis.line = element_line(color='black'),
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank())+geom_vline(xintercept =zstar, 
                color = "red", size=0.3)

dev.off()





##replication of panel b##
a=14000;
b=51500;
zstar=14000;
x0=500;
binwidth=x0;


y0=5;



bins_excl_r <-1
         
            firstpass_prep <- bunching::prep_data_for_fit(binned_data0813,zstar =zstar,  binwidth = x0,
                          bins_l =a/x0, bins_r=b/x0,rn=c(5000,10000), bins_excl_r=bins_excl_r,bins_excl_l=0/x0,poly=y0)
            firstpass <- bunching::fit_bunching(firstpass_prep$data_binned, 
                firstpass_prep$model_formula, binwidth=x0, notch=TRUE)
            B_below <- firstpass$B_zl_zstar
            M_above <- -firstpass$B_zstar_zu
            if (M_above > B_below) {
                stop("Missing mass above zstar is larger than bunching mass below. Are you sure this is a notch?")
            }
            available_bins <- (max(firstpass_prep$data_binned$bin) - 
                zstar)/binwidth
            DoF_remaining <- nrow(firstpass_prep$data_binned) - 
                nrow(firstpass$coefficients) - 1
            notch_iterations_bound <- min(available_bins, DoF_remaining)
            zu_bin <- bins_excl_r
            tmp_firstpass_prep <- firstpass_prep
            while ((B_below > M_above) & (zu_bin < notch_iterations_bound)) {
                zu_bin_final <- zu_bin
                zu_bin <- zu_bin + 1
                newvar <- paste0("bin_excl_r_", zu_bin)
                tmp_firstpass_prep$data_binned[[newvar]] <- ifelse(tmp_firstpass_prep$data_binned$z_rel == 
                  zu_bin, 1, 0)
                   tmp_firstpass_prep$model_formula <- stats::as.formula(paste(Reduce(paste, 
                  deparse(tmp_firstpass_prep$model_formula)), 
                  newvar, sep = " + "))
                tmp_firstpass <- bunching::fit_bunching(tmp_firstpass_prep$data_binned, 
                  tmp_firstpass_prep$model_formula, binwidth=x0, notch=TRUE)
                B_below <- tmp_firstpass$B_zl_zstar
                M_above <- -tmp_firstpass$B_zstar_zu
                if (B_below > M_above) {
                  firstpass_prep <- tmp_firstpass_prep
                  firstpass <- tmp_firstpass
                  zu_bin_final <- zu_bin
                }
            }

mb<-zu_bin*x0+zstar;  
plot_binned_data0813<-tmp_firstpass_prep$data_binned;

plot_binned_data0813$cf<-tmp_firstpass$cf_density;

plot_binned_data0813$pfreq<-100*plot_binned_data0813$freq/sum(plot_binned_data0813$freq);
plot_binned_data0813$pcf<-100*plot_binned_data0813$cf/sum(plot_binned_data0813$freq);



pdf(width = 8, # The width of the plot in inches
    height = 4) # The height of the plot in inches


ggplot(plot_binned_data0813)+geom_line(aes(x=bin,y=pfreq),color="black",size=0.4)+geom_line(aes(x=bin,y=pcf),linetype="longdash",color="red",size=0.4)+ 
  labs(x="Loan Amount",y="Percentage of Loans")+geom_vline(xintercept = 13880, linetype="dotted", 
                color = "black", size=0.3)+
scale_x_continuous(breaks=c(5000,15000,25000,35000,45000,55000,65000))+
geom_text(parse=TRUE,size=3,aes(x=13300, label ="italic(K[L])", y=8), colour="black", angle=0, text=element_text(size=0.3))+geom_vline(xintercept = mb, linetype="dotted", 
                color = "black", size=0.3)+geom_text(parse=TRUE,size=3,aes(x=21300, label ="italic(K[U])", y=8), colour="black", angle=0, text=element_text(size=0.3))+ theme_bw() +
  theme(axis.line = element_line(color='black'),
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank())+geom_vline(xintercept = zstar, 
                color = "red", size=0.3)

dev.off()



##replication of panel c##

a=25000;
b=40500;
c=0;
zstar=25000;
x0=500;
binwidth=x0;


y0=5;

bins_excl_r <- 1
         
            firstpass_prep <- bunching::prep_data_for_fit(binned_data14200,zstar = zstar, binwidth = x0,
                          bins_l =a/x0, bins_r = b/x0, rn=c(5000,10000), bins_excl_l=c/x0,bins_excl_r=bins_excl_r,poly=y0)
            firstpass <- bunching::fit_bunching(firstpass_prep$data_binned, 
                firstpass_prep$model_formula,binwidth=x0,notch=TRUE)
            B_below <- firstpass$B_zl_zstar
            M_above <- -firstpass$B_zstar_zu
                       if (M_above > B_below) {
                stop("Missing mass above zstar is larger than bunching mass below. Are you sure this is a notch?")
            }
            available_bins <- (max(firstpass_prep$data_binned$bin) - 
                zstar)/binwidth
            DoF_remaining <- nrow(firstpass_prep$data_binned) - 
                nrow(firstpass$coefficients) - 1
            notch_iterations_bound <- min(available_bins, DoF_remaining)
            zu_bin <- bins_excl_r
            tmp_firstpass_prep <- firstpass_prep
            while ((B_below > M_above) & (zu_bin < notch_iterations_bound)) {
                zu_bin_final <- zu_bin
                zu_bin <- zu_bin + 1
                newvar <- paste0("bin_excl_r_", zu_bin)
                tmp_firstpass_prep$data_binned[[newvar]] <- ifelse(tmp_firstpass_prep$data_binned$z_rel == 
                  zu_bin, 1, 0)
      tmp_firstpass_prep$model_formula <- stats::as.formula(paste(Reduce(paste, 
                  deparse(tmp_firstpass_prep$model_formula)), 
                  newvar, sep = " + "))

   tmp_firstpass <- bunching::fit_bunching(tmp_firstpass_prep$data_binned, 
                  tmp_firstpass_prep$model_formula,binwidth=x0, notch=TRUE)

                B_below <- tmp_firstpass$B_zl_zstar
                M_above <- -tmp_firstpass$B_zstar_zu
                if (B_below > M_above) {
                  firstpass_prep <- tmp_firstpass_prep
                  firstpass <- tmp_firstpass
                  zu_bin_final <- zu_bin
                }
            }
     
mb<-zu_bin*x0+zstar

plot_binned_data1420<-tmp_firstpass_prep$data_binned;

plot_binned_data1420$cf<-tmp_firstpass$cf_density;

plot_binned_data1420$pfreq<-100*plot_binned_data1420$freq/sum(plot_binned_data1420$freq);
plot_binned_data1420$pcf<-100*plot_binned_data1420$cf/sum(plot_binned_data1420$freq);


pdf(width = 8, # The width of the plot in inches
    height = 4) # The height of the plot in inches


ggplot(plot_binned_data1420)+geom_line(aes(x=bin,y=pfreq),color="black",size=0.4)+geom_line(aes(x=bin,y=pcf),linetype="longdash",color="red",size=0.4)+ 
  labs(x="Loan Amount",y="Percentage of Loans")+geom_vline(xintercept = 24880, linetype="dotted", 
                color = "black", size=0.3)+
scale_x_continuous(breaks=c(5000,15000,25000,35000,45000,55000,65000))+
geom_text(parse=TRUE,size=3,aes(x=24200, label ="italic(K[L])", y=11), colour="black", angle=0, text=element_text(size=0.3))+geom_vline(xintercept = mb, linetype="dotted", 
                color = "black", size=0.3)+geom_text(parse=TRUE,size=3,aes(x=42200,label ="italic(K[U])", y=11), colour="black", angle=0, text=element_text(size=0.3))+ theme_bw() +
  theme(axis.line = element_line(color='black'),
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank())+geom_vline(xintercept = zstar, 
                color = "red", size=0.3)

dev.off()
